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Abstract 

In this paper we derive time reversal imaging functionals for two strongly causal acous- 
tic attenuation models, which have been proposed recently. The time reversal techniques 
are based on recently proposed ideas of Ammari et al for the thermo-viscous wave equa- 
tion. Here and there an asymptotic analysis provides reconstruction functionals from first 
order corrections for the attenuating effect. In addition, we present a novel approach for 
higher order corrections. 

1 Introduction 

Photoacoustic Imaging is a promising imaging technique for visualizing biological material 
parameters. In experiments, the medium is exposed to a short pulse of an electromagnetic 
wave. The medium absorbs a fraction of the induced energy, heats up, and reacts with 
thermoelastic expansion. This in turn produces acoustic waves, which can be recorded and 
which are used to determine the electromagnetic absorption coefficient |18|. El \TT\ [19] . These 
coupling properties explain why photoacoustic imaging is referred to as hybrid or Coupled 
Physics Imaging. For some recent progress in hybrid imaging we refer to the surveys [14] [4] . 

In this paper we investigate the method of time reversal in attenuating media as it was 
introduced in [21 \T7\ [3] for the thermo-viscous wave equation. In these references, the goal is 
to construct a parameter dependent family of approximate reconstruction functionals, which 
allow for approximation of the initial datum of the thermo-viscous wave equation. For time- 
reversal technique for imaging in general, see for example [HEIH]. 

In this work, we are using the asymptotic techniques from [21 [T7] to develop time reversal 
algorithms for the Nachman-Smith-Waag (NSW) |15j and the Kowar-Scherzer-Bonnefond 
(KSB) models [12], [13] . These two models satisfy a strong causality property. That is, the 
solutions of these equations are zero before the initialization and satisfy a finite front wave 
propagation speed. We emphasize that the thermo-viscous wave equation considered in [21 
[17] is not strongly causal [9]. While there is a large literature on inversion formulas and 
time reversal algorithms for the standard wave equation much less is known in the case of 
attenuating media [5] [12] [2] [15] . Partially, this is due to the fact that, so far, the reference 
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attenuation model (that is, the governing wave equation) has not been established in this 
field. 



2 Review on time-reversal for the thermo-viscous equation 

Let Q be a bounded domain in R 3 . Standard photoacoustic imaging consists in determining 
the absorption density /, which is assumed to have compact support K in Q, in the acoustic 
wave equation 

§0M) - A P (x,t) = ^(t)f(x), (x,t) e R 3 x 1, 

9t dt (2.1) 

p(x,t) = and ^(x,t) = Q, x G M 3 , t < 0, 

from measurement data g(y, t) := p(y, t) for some y G d£l and t G [0, T], where T is supposed 
to be sufficiently large. 

In the following, let a := a(x) be a positive function, describing attenuation. The imaging 
method considered in [21 [17] consist in reconstructing / from data 

g a (y, t) := p a (y, t) for all y G dQ, and i G [0, +oo). 

where p a solves the thermo-viscous wave equation, 



I ( M )_ (x + a -) Ap a (x,t) = -±(t)f(x), (i,t)£l 3 x 



dt 2 v ' ; V 94 7 ' 9t 



p a (a;,f) = and ^-(x,i) = 0, i£R 3 , i«0 



(2.2) 



To derive the imaging technique we use rj(x, •), the fundamental solution at x of the 
Helmholtz equation, 

oj 2 u(x, y) + (1 + iao;)Aj / n(a;, y) = -<5*(|/), y G K 3 . (2.3) 

In [21 [T7] the following results have been shown: 

• For fixed cut off parameter p > the function 

v a s , P (x,t) := f [ iojT a UJ (x,y)g a (y,T - s)da(y)e-^ t - s Uu , (2.4) 
2vr J_ p J dQ 

satisfies the thermo-viscous wave equation 

9 v ( d \ dS 

-r^ 9a(x,T - s)S 9 a, (2.5) 



Op(x, t)-\I-a—) Av(x, t) = S p 



where S p is as in (|5.21|) . In this formula the necessity of the regularization becomes 
evident because in the unregularized form the right hand side consists of a product of 
two distributions, which a-priori is not well-defined. 
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• Moreover, for small a it follows from the results in [17] that 

l a p (x)= [ T v a (x,T)ds — ► f(x), asp^+oc. (2.6) 
Jo 

Moreover, [171 Remark 2.3.6] also gives a reference how to obtain appropriate values 
for the cut off parameter p. The regularization is required so that the function v% p (x, t) 
is well-defined in (I2.4p . 



3 The KSB-model 



In this section, we are deriving an imaging functional for the KSB model [12, Eqs.(28),(88)] 
following the time reversal approach of [2] for the thermo-viscous wave equation, as outlined 
above. 

Let ao > and 7 E (1,2] be fixed. Then the KSB model assumes that the attenuated 
pressure p a satisfies the equation 



(a Z + L 1/2 



2 0V/ X t A a / \ T 9S (t) Pl s 
% t)-L Ap a (x, t)=L 



dt 2 



dt 



x E 



t E 



dp a 



p a (x,t) = and ^-(x,t) = 0, x E M 3 , t < 0, 



1 



(3.1) 



-l\l/2 



where L 1 / 2 is the convolution operator (in time) with kernel T 1 (l + (— iro^) 7 

v27r L 

- we emphasize that L = L 1 / 2 o L 1 / 2 . Here Dj~ l denotes a fractional time derivative operator 



of order 7 — 1 PjQJ 

The Fourier transforms, p a := F[p a ](oj) and p := T\p\{ui), of (|3.ip and (|2.ip satisfy the 
Helmholtz equations: 

K 2 {oj)p a {x) + Ap a (x) = icjf(x), and oj 2 p(x) + Ap(x) = ±uif(x) , 



respectively. Here 



k(uj) = OJ I 1 + 



x (l + C-i^-ijvay • 

Using the particular form of the Helmholtz equations it follows that 

OJ 



(3.2) 



V 



k(oj) 



7"[p](kH), 



which after applying the inverse Fourier transform yields 

p a (x,t)=C a [p(x,-)](t) 

with 



0J 



k(oj) 



1 



(*) = sz 



a; 



2vr 7 R k(w) 



-iwt /" e^M-^^rfady . (3.3) 
Jm 
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The goal is to find an imaging operator of the form, that is functions R and A, 

a; A (a; 



2ir J R k{oj) J R 



(t), (3.4) 



k{uj) 

for which the following expansion with respect to «o ( as in (|3.2p ) holds: 

C* a C a [4>]{t) = 4>(t) + o(a ), asao^O. (3.5) 

Thereby, 

Z*a[<l>]{t) ■= ^ I ^rA( W )e iS ^ / e-^^dsdu (3.6) 

is the adjoint of C a . 

The derivation of imaging operators as in (|3.4p follows general principles, which are used 
later on for the other attenuation models as well. The principle consists in construction 
auxiliary functions Ai, ui, and U2 which are determined from k by general construction 
without taking into account the special structure of the function. The explicit construction 
comes at a later stage. We introduce the auxiliary function Ai : R — > R which is related to k 
in the following way: 

k(u) =w(l-aoAi(w)) + 0(ag). (3.7) 
Then, it follows by expansion with respect to ao that 

to . . , , „, 2 



1 + aoAi(u>) + O(a ) for fixed uj G 



k(oj) 

e iK^)s = e ius^ _ iQQa;Al ( w ) s ) + 0(^2) for fixed g £ R ^ 

and consequently for fixed s,oo G R 

!«(«). = e i^ s (1 + ao ( Al ( w ) _ iujX^s)) + 0(ag) . (3.8) 



-e* 



k(w) 
This shows that 

£«[0(a)](t) = 0(0 + ao^" 1 [Ai(w)^[0(*)](o;) - iwAi(o;)J*[*0(s)](a;)] (t) + o(a ) 

= £ o [0(s)](t) + "o^ 1 [Ai(w)^[^(«)](a;) - iwAi(w).F[^(*)]M] (t) + o(a ) . 

Note that for a = 0, Cq = X. In order to get an explicit form of C a we introduce auxiliary 
functions Vi : R — >■ R, i = 1, 2, and 



re(o;) := w(l - a ^i(w)) + , 
A(w) := 1 + a ^(^) + O(ao). 



(3.9) 



1 In [2] [17], (k(w) — w) is an imaginary function, therefore the Stationary Phase Method is used for the 
asymptotic analysis of the previous integral operators. Here the relevant function, ( — aoAi(o;)) in its general 
form, is complex (with non- vanishing real part) and the usage of the method of Steepest Descent would provide 
the relevant asymptotic expansion. The Stationary Phase Method can be seen as a special case of the method 
of Steepest Descent. In this work, the relevant real part is zero for the NSW and the thermo- viscous model and 
consequently the Stationary Phase Method can be applied. However for the KSB model this part is non-zero 
and the use of the method of Steepest Descent is needed. 
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Again, by expansion with respect to ao it follows that 

1 + a (^(w) + 1/2M) + O(«o) . 



K{OJ) 

and consequently, 



K(w) 



(3.10) 



Therefore, we have 



+ ao [(n(-w) + M-u))T[cf>(s)}(u)] (t) 
+ o(a ). 

The goal is to determine 1/1,1*2 such that the corresponding operators C a and C* a satisfy 
(13. 5h . Using the two expansions for £ a and £* it follows that 

£:>c a ^( 5 )](t) =<t>(t) 

+ a j^ 1 [(Ai(w) + i^(-w) + ^(-a;)) ^[0(s)](w)] (t) 

+ [-iwAjCwJ^WKw)] (t) (3.11) 
+tJF- 1 (*)} 

+ o(a ) 

To satisfy (|3.5p we require that the first order term in ao of the equation (|3.1ip has to vanish. 
By taking the Fourier transform of this term we see that the term vanishes if 

(AiM + uxi-oj) + v 2 (-u)))T[4>](u) - wAiCw)^^^) + A. (uj Vi {-u)F[4>]{uj)) = 0, (3.12) 

where we have used the property (|5.22p with n=l. Now, it is straightforward to see that the 
solution of the following system 

M „ )+ ^ )+ ^ )+ fcfc^_„, (3i3) 

— uj\\{lS) + ujvi{— uj) = , 

satisfies the equation (|3.12p which directly implies (|3.5p . Equivalently, we get the following 
conditions 

v\(— uj) = X\(oj) and ^(—oj) = — 3Ai(w) — u — y- — - . (3-14) 

duo 

Now, we introduce T w (x,y) and T^(x,y) which are the fundamental solutions of the 
Helmholtz equations 

uj 2 Tjx,y) + A y T ul (x,y) = -5 x (y), y £R 3 (3.15) 
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and 

k(u;) 2 r a ul (x,y) + A y r a ul (x,y) = -\(L J )5 x (y), 
respectively. We can prove that 



y e 



dT 

M 



where C a is defined in (|3.4p . 



F(x,y,t,r) = T 1 {r w (x,y)} (t - r) 



and 



f a (x )2 /,t,r) =^ 1 {C(x,y)}(t-r), 

and J 7-1 denotes the inverse Fourier transform with respect to oj. 
Then, we define the function Vg(x,t) by 



(3.16) 
(3.17) 

(3.18) 
(3.19) 



vt(x,t) = -±- [ [ iut a UJ (x,y)ga(y,T-s)da(y)e- i ^ t -^ ( Lj, (3.20) 
£k Jr Jan 



where we recall that g a (y,t) := p a (y,t) for all y G dQ and t G [0, T). 
For the KSB model, it follows from fl32J) and flO) that 



A 1 (o;) = -(l + (-ir ^- 1 ) 



Using this in ()3.14j) we get 



i/iM = - (1 + (iro^ 1 ) 



1N-1/2 



and 



7 ~7 /-, , ^_ . w-n-l/2 ,7-1 



1 + (ir w)' 



+ 



'1 + (ir w) 



7-l\-3/ 2 



2 v v " ' ' 2 
Using these expressions for v\{uS) and ^(w) we suggest the following choice for k(cj) and 
A(w): 

R{uj) = oj(1 — a.QV\{uSj) and \{oj) = 1 + ao^2(w) , (3.21) 

respectively. Note, that there is no remainder term in ()3.14p . 

By applying the expressions (|3.2ip into (|3.16p , then the previous definition of the function 
Vg(x,t) yields the following identity: 



L 1 ' 2 L + ao((7-l)X+(7- 7 )L) -£ (g a (x,T - s)d dQ ) forxGfi. 



06, 



Here 
and 



dt 

L=i+(- T0 y- 1 Dr\ 

g a (y, t) := p a (y, t) for all y G dfl , t G [0, T] , 



(3.22) 
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where T is supposed to be sufficiently large such that p a (x, t) = = ^-(x, t) for t > T and 
x G Q. 

In the following subsection we prove that the functional 



l a (x) 



T 



v a s {x,T)ds 



is an approximation of the initial datum f(x). For doing this we have to make some regular- 
ization of the relevant operators. 



3.1 The regularized time reversal functional 

Matters of convergence of some infinite integrals defined above, suggest a regularization in 
the same way as in [2] and [17]. We define the function 



< P (x,t) 



1 /•-" 
2^ 



-p J an 



iuf Ux,y)g a (y,T- s)da(y)e- i ^ t - s Uu;. 



(3.23) 



as an approximation of Vg(x,t) defined in (|3.20p . Moreover, we define the regularized funda- 
mental solution of r a as in (I3.19h 



f a Jx,y,s,t) 



2tt 



the regularized operator C ap defined in (|3.3 



£aJ4>]{t) 



1 

2^ 



;A(o 



-p k(u) 



and its adjoint 



27r ./_„ R{oj) 



4>(s)dsdui. 



- P <«\>~, Jo 

Using these definitions we write the approximated version of the equations (|3.5p and (|3. 1T[) , 
respectively, that is 



C* Calm) = Sp[m)+o(a ), 



and 



£ — r 

dt a " 



dr 



dt 



(3.24) 



(3.25) 



where the operator S p [(f)] and the function T were defined in (|5.2ip and (|3.18p . respectively. 

Similarly to the previous subsection, applying the definition (|3,23p of (x,t) in equation 
(|3.16|) with usage of the expressions (|3.21|) and (|3.21|) we obtain the following wave equation 



(a X + L^) 2 ^ - ft* L a) t) 
■S P [ L + a (( 7 - 1)1 + (7 - 7)1) ) ^ 



(3.26) 



(g a (x,T - s)5 dn ) for x e VL . 
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Finally, we can obtain the reconstruction functional I™. Indeed, since equations (|3,24|) and 
(|3.25p hold, then Proposition 2.3.5 in [IT] shows that 

Z?(aO= / v s,p( x iT)ds — > f(x), asp ^+00. 
J 

Remark 1. The latter proposition suggest that the larger p we choose, the better approxi- 
mation we get. However, the previous computation of Tp(x) involves the integration of the 
fundamental solution T^(x,y) which grows exponentially as exp {9{k(u;)}|x — y\}. In order 
to ensure stability of X^(x) this term must not be greater than one \1T[ Remark 2.3.6]. For 

large u, the expression \$${k(u))}\ is behaving like (and is less than) cxo\u)\ sin ^ ^ , 7 6 (1, 2]. 
So one should not use frequencies larqer than — ,. 1 . Hence, we qet p ~ — ,. 1 , n . to be 
the threshold for the imaging functional stability, where diam{£l) denotes the diameter of 
the domain fi. A finer estimation of the threshold can be given if we use that is 



behaving like ao T o 2 M 2 7 sin ^ 7 ^ w , 7 £ (1,2], for large values ofui. Consequently, 



7-1 

T 3-7 



2 

a,odiam(Q) sin ( 7 ~ 1 ) 7r 1 3 7 



■1 



twtft 7 e (1,2]. 

4 The NSW model 

Let p a satisfy the following problem 

d \ d 2 p a / 9\ 4 „ / 3 \ d<5 , 

X + T atJ ( X + T ^J AP = ( X + T ^J ^ (41) 

along with the conditions in ()3.ip . where we consider the NSW model for one relaxation 
process, as defined in [12]. Moreover, we assume that r > f > so that the strong causality 
condition in |12j is satisfied and that r and f are small and of the same magnitude. 
Substituting f = a^f and r = a^r, we find 



, , /l — iwf 1 — iaour 

«M = W V 1 : =UJ \ S : • 

V 1 — iwt VI - iao^r 

By applying the expansion (13. 7j) for /t(u;) in terms of ao we get 

r — r 

Ai(w) = -iw— . 

Consequently, the auxiliary functions i/j, i = 1, 2 are given by the conditions (|3.14p and read 
as follows 

v\{uj) = iu — - — and V2{u) = — 2iw(r — f) . (4.2) 
Therefore, we obtain the expansion of k(cu) and X(to) from (|3.9p . i.e., 

k(oj) := u(l - a vi(u)) + 0(a%) , 
A(w) := l + Q ^ 2 (w) + C , (ao), 



S 



where now the functions vi , i = 1 , 2 are given by (|4.2p . 

We make the following choices for the functions k(ui) and A(w) 



2 



/l + icjf / 1 + iwr . 
k(u) = uj\ ; and \{uj) = I : , (4-3) 



1 + icoT ' \1 + iwr 

which satisfy (13. 9p for non-vanishing O(ag) terms. From here, using the Helmholtz equa- 
tion (|3.16p and with the same arguments as in the previous section we derive the following 
regularized time reverted attenuated equation (which is of course not unique) 




-Sp 



T dt at 



a (x,T - s)5 9 n) , 



(4.4) 



where S p is defined in (15.2ip , Here 

g a (y, t) := p a (y, t) for all y G d£l , t G [0, T] , 

where p a (y,t) satisfies equation (|4.ip and T is supposed to be sufficiently large such that 
p a (x,t) = = %{x,t) for t > T and iGO. 

The reconstruction imaging functional is given by (|2.6p . i.e. 

Zp( x ) = / v s P i x iT)ds — > f(x), as/5-^+oo. 

Remark 2. For i/ie case o/ N relaxation processes, the procedure will be conceptually the 
same. Now, the attenuated wave equation has a more complicated form fTblj and we have 



i N t 

1 v-^ 1 — lWTi 



N ^ 1 - iuT-j 
j=l o 



Here, we assume that {Tj,fj}i are small and of the same magnitude, i.e. all {Tj,fj}^ are 
of order 0(ao). Then the expansions of the functions k(oj), k(uj) and A(w) in terms of oq, 
as they were given in (|3,7p and (|3,9p , and the relevant asymptotic analysis allow us to make 
(as previously) the following choices 



i 



1 N 1 , • 

1 1 + 1WT. 



N 




which lead to the corresponding time-reverted attenuated equation. One can find the corre- 
sponding wave equation by applying the inverse Fourier transform on the Helmholtz equation 
(|3.16p . for the later values of k{uj) and X(uj). This procedure will lead to an equation similar 
to (|4.4p . but this time one obtains a more complicated form. 

Following the analysis of the previous section we will choose the value of the truncation 
parameter p appearing in (14. 4p by finding the behaviour of the expression \^s{k{ui))\, with k(uj) 
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UJ 2 

given in (14, 3p . For \to\f < \uj\t < 1 we find that \*Q(k(lj))\ is behaving like — (r — f), which 
(following the arguments from Remark 2.3.6]) yields the truncation parameter 

1 

P 



\/ diam(£l) (r — f) 

Moreover, this value of the truncation parameter is sufficient for the other two cases, i.e. 
for 1 < \uj\t < \uj\t and \oj\t < 1 < \uj\t, when \^s(k(lj))\ is behaving like — 7) arid 




,,_ - 2 y.-r ,,\r 1 ' respectively. 

For the case of N relaxation processes, one can use similar arguments to get an estimation 
of the truncation parameter p. In the case that \u)\fj < \ui\rj < 1, for all j = 1, . . . , N we find 

UJ 2 

that \Q(k(uj)}\, with k{u) given in (|4.5p . is behaving like —— Ylj=i( T j which yields the 

truncation parameter 



N 



p ~ 



For the several other cases, the arguments of the previous remark, along with the usage of the 
triangular inequality, give us the opportunity to observe that the above-mentioned estimation 
for the truncation parameter is sufficient. 

Remark 3. The formal procedure outlined above also applies to thermo-viscous model, as 
this was defined by the wave equation (84) in fT2^ . This is the special case of the NSW-model 
with one relaxation process and f = 0. Note that the thermo-viscous wave equation (84) in 
jll2f refers to a not strongly causal model and has a different RHS from the thermo-viscous 
equation (12. 2p . which was analysed in \1 7| /. 



5 Higher order terms 

In this section we describe the procedure for evaluating higher order terms of the operators 
C a and £*, defined in (]3.3h and (I3.4D and consequently a higher order approximation of the 
reconstruction functional X". 

In particular, this method allows determining the higher order terms of the asymptotic 
expansion for the functions and A(w), appearing in (|3.4p . We make the following ansatz 

00 

k(u) = uj \j(uj)a j , Ao(w) = 1 (5.1) 

and consequently 

00 

^ = £m;("K, (5-2) 

with 

Po(uj) = 1, MiM = fi2(cj) = - A 2 (w) (5.3) 
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and in general 

lij(u) = Vj ({Afc(w)}Li 

where Vj denotes a polynomial (on several variables) of order j. In addition, the expansion 
(|5.ip yields the following 



j=0 

with 

1 2 

tp (uj,s) = l, ipi(oJ, s) = —iuj\i(u)s, ip2(u, s) = iwA 2 (o;)s + - (iwAi(a;)s) (5.5) 
and in general 

Y>j(w,s) = Qj ({iu}\ k (uj)s}{ =1 ^ , (5.6) 

where Qj denotes a polynomial (on several variables) of order j. In the case of higher order 
terms we make the ansatz 

oo 

«(w) = u ^2(-l) j Xj{-uj)aK (5.7) 
j=o 

This expansion is without loss of generality. In the procedure described in section [3l concern- 
ing the first order approximation of k(lj), the relevant terms were considered unknown and 
one had to determine them. However, applying the previous expansion we state a consistency 
condition for the time-reversal algorithm which, instead of giving a system of equations, yields 
a set of identities, as we will see below. 

In the same way with the previously mentioned expansions we get 

oo 

OJ 



(5.8) 

and 



e i*( w > = e ^'J2i>3(u,8)aP, (5.9) 
j=o 

with 

^(u, s) = Qj ({iwA fc (-w)s}{ = i) , (5.10) 



where Qj denote the same polynomials as in (|5,6 
Finally, we consider the expansions 



ah = ]T/?>k, A)M = i (5-11) 

3=0 



and 



A(w) = ojj^^ai , 70 (o;) = 1. (5.12) 



V ^ 3=0 
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These, expansions along with (15. 7|) yield the following relation 



(5.13) 



Now, since we know explicitly our strategy consists of determining A(w); conse- 

quently the Helmholtz equation (|3.16|) . under the procedure described in the relevant sec- 
tion, provides the corresponding time-reverted wave equation. So, in what follows out target 
is to describe a procedure to find f3 n (u), n G N, equivalently determine the terms of the 
asymptotic expansion of A(u>). 

Applying the previous expansions in the expressions (13.3H and (13. 61) . we get 



(5.14) 



£a [</>](*) = and C*M(t) = 



fc=0 



fc=0 



with /o = go = Id, {/fc> <?jfc}fc = i being operators that can be obtained explicitly in terms of 
{Xj, ldj,ipj,tpj,^j}j =1 , for fc = 1, . . . , n, respectively. In particular, {/&,fffc}fe=i are obtained 
explicitly with use of the expressions (|5.3p . (|5,5p . (|5.10p . i.e., 

A [<£(*)](*) = [-iwAiCw^^WKw) + Ai(w)^[0( S )](a;)] (t), 



/ 2 [^( S )](t)=J- 1 



;A?( W ) - A 2 M) W)K") - iw (Af (w) - A 2 (w)) ^(a)](o;) 



+-(ia;A 1 (a;)) 2 ^ 2 ^( S )](a;) 



(*), 



S2[^(fl)](t) = .F" 1 [ 72 (-o;W( S )]H] [io;(7 1 (-a;)A 1 (a;) - A 2 (w))^(a)](a;)] (t) 



+i 2 F" 1 



-(ia;Ai(a;)) 2 ^[^(a)](a;) 



(t). 



Consequently, we get the following asymptotic expansion 



a 2 +o(a 2 ), (5.15) 



when a — > 0, with {/&, <7fcH=i defined above. The identity 

CC a [0](i) = 0(t)+o(a n ), a^O (5.16) 
is satisfied, up to the 1st order (equivalently the identity (|3.5p is satisfied), when 

T[fM+gi[(t>}}(oj)=Q, V W £K. 



Using the definitions of the operators fx and g\ given above and the property (15.22P this 
condition yields the following equation 



7i(— u) = — 2Ai(w) — oj 



d\i(u}) 



(5.17) 
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which is equivalent with the conditions (13. 13f) . By applying the condition (15.17P in the 
expression (|5.13p for n = 1, we get that 

= -x^u) + 71 (- W ) = -3Ai(w) - ^^J^, (5-18) 

which is the first order approximation, of the function \(cj), as this was defined in (|5.1ip . 
Prom the expressions (|3.9p and (|5.1ip . one can see that V2 = Pi, as they were introduced 
in these equalities, respectively. So, equation (|5.18p is exactly the second of the conditions 
(I3.14p . Note that the first of these conditions, is satisfied as an identity by the choice of the 
expansion of k(u), given in (|5.7p . 

Following the same procedure for the 2nd order in the identity fj5.16[) . the expression 
(I5A51) yields 

^[/ 2 [0]+3 2 [0]+3l[/i[0]]]M = O, 

which after some calculations where we made use of the definitions of the operators /2, #2 
and the fact that 



-1 



iwA5(w)^(s)](w) - (iwAi(w)) 2 ^(s)](w) (t) 



we obtain the following equation 

72 (-w) = -A?(w) + A 2 (w) - 7i(-w)Ai(£j) 

d / , >9 ,\ 1 d 2 /, , , nn2 \ (5-19) 



^ du) 



L (-A?(w) + A 2 (w) - 7i(-w)Ai(w)) ) - ^ (MiH) 2 ) 



Substitution of the conditions (|5.17p and (|5.19p in the expression (|5.13p provides the second 
order approximation for the function A(w), i.e. /3 2 (w), as this was defined in (|5.1ip . 

A general treatment of this problem would appear as follows: The general form of (|5.15p 
would be 

n 

C* a C a [<P](t) = <t>{t) + Y, h k [<t>]{t)a k + o(a n ), a -> 0, 

k=l 

where {hk}^ =1 are operators expressed explicitly (after some calculations and the use of the 
previous expansions) in terms of {Ay, 7y}| = i, for k = 1, . . . , n, respectively. Consequently, the 
condition 

T[h k [<t>]]{u) = 0, Vw£R 

yields 7 fc(— u) as an explicit expression of {Aj(o;)} - =1 , { 7 j(— w)}^~| and their derivatives up 
to order fc. Hence, the expression f)5 . 13|) yields the fc-th order approximation for the function 
A(w), i.e. h(u)). 
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Appendix 

In this paper we use the basic notation: 

• 5 S denotes the one-dimensional Dirac distribution at s € ML Moreover, we denote by 6 X 
the 3-dimensional (5-distribution with center x - x is always a 3D variable. 

• T denotes the Fourier transform 



• If we want to specify the argument, with which respect the Fourier transform is applied 
it is specified as follows: 




(5.20) 




In this case it is with respect to the variable s. 




(5.21) 



• 5qq, denotes the surface Dirac mass on 8ft. 



Moreover we use the following property of the Fourier transform 



n 



du n 



(5.22) 



for n 6 N. 
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